% To compute the latitudinal profile of B_x (ionospheric) from CM4 for
% all the selected mjd instances of satellite and observatory pairs at all
% the equatorial observatories. Fit a polynomial and get the values at each
% station pai
%latest date Feb 10 2006

%observatory locations GGLAT GGLONG

ABG	=[18.63	72.87]; %1
ETT	=[9	78];        %2
HYB	=[17.42	78.55]; %3
TIR =[8.67 77.82];  %4
HUA	=[-12.05 284.67];%5
FUQ	=[5.47 286.27]; %6
MBO	=[14.39 343.04];%7
GUI	=[28.32	343.56];%8
AAE= [9.02 38.77];  %9
QSB	=[33.87	35.64]; %10
GUA	=[13.58	144.87];%11
CBI	=[27.1	142.18];%12
PND =[11.92	79.92];
ELT = 
YAP = 


%selected observatories

%ETT HUA MBO AAE GUA

location = [9 78;-12.05 284.67;14.39 343.04;9.02 38.77;13.58 144.87];

%load the dH files from observatories


%load index files

load('c:\manoj\projects\eej\Sq-India');
phi = 79;


%generate the variable vector for firing the CM4
load c:\manoj\projects\eej\mjdkpf107 mjd kp F107;
path(path,'c:\manoj\geomag\mf\cm4\');
r = 6371.2;
Dst = -50;
options.n_min = [ 1 14]; % n_min
options.n_max = [13 45]; % n_max
options.cord = 1;     % geocentric coordinates
options.pred = [1 1 1 1 0 0];
options.filename =  'c:\manoj\geomag\mf\cm4\umdl.CM4';

theta = [40:2:140];

pp = 1
for i = 1:length(theta),
    for j = 1:length(index),
      mjd1(pp) = mjd(index(j));
      theta1(pp) = theta(i);
      phi1(pp) = phi;
      Dst1(pp) = Dst;
      F1071(pp) = F107(index(j));
      r1(pp) = r;
      pp = pp+1;
    end;
end;



[B_x1, B_y1, B_z1] = cm4(mjd1, r1, theta1, phi, Dst1, F1071, options);

%Polynomial fitting of the Sq B_x fields 
load('C:\Manoj\projects\eej\OBSdH1');
eejdelta = reshape(squeeze(OBSdH(2,:,:))',[1,1096*24]);
eeodelta = reshape(squeeze(OBSdH(3,:,:))',[1,1096*24]);

EEJ = ETT;
EEO = HYB;

theta = [40:2:140];

B_x = reshape(B_x1,[length(index) length(theta) 7]);
B_y = reshape(B_y1,[length(index) length(theta) 7]);
proj = [B_x(1,1) B_y(1,1)]/(sqrt(B_x(1,1)^2+ B_y(1,1)^2));

for i = 1:length(index),
B_h = (B_x(i,:,5)+B_x(i,:,6))*proj(1) + (B_y(i,:,5)+B_y(i,:,6))*proj(2);
% B_x1 = (B_x(i,:,5)+B_x(i,:,6));
[p,s] = polyfit([40:2:70 100:2:140],B_h([1:16 31:51]),5); 
%[p,s] = polyfit([40:90 120:140],B_h([1:50 80:101]),8); %HUA - FUQ
%sq_est(i,:) = polyval(p,[90-EEJ(1),90-EEO(1)]);
%
y = polyval(p,theta);
% y1(i,:) = y;
% B_x2(i,:) = B_x1;
plot(theta,B_h,'r.-');
% plot(90-[40:1:140]-8,B_h);
% set(gca,'XDir','reverse');
hold on;
plot(theta,y,'k');
% plot(a(:,1),a(:,2),'k');
plot(90-EEJ(1),eejdelta(index(i)),'r*',90-EEO(1),eeodelta(index(i)),'r*');
plot(90-EEJ(1),eejdelta(index(i)),'r*',90-EEO(1),eeodelta(index(i)),'r*');
plot(90-EEJ(1),eeodelta(index(i))+(sq_est(i,1)-sq_est(i,2)),'k*');
 [iy, im, id, ih, min, sec] = datevec(mjd(index(i))+730486); % t is MJD2000
  title(sprintf('ETT-HYB %d /%d /%d, %d:%d UT',iy,im,id,ih,im));
% 
 pause;
 hold off;
%
end;

%{
%load solarflux data
[a11,b11] = xlsread('c:\manoj\geomag\indices\solarflux.xls');

% kp 

load c:\manoj\geomag\indices\kp_XXXX.list;
global data_ap;
data_ap = kp_XXXX;
clear kp_XXXX;


% mjd, F107

st1 = date2jd(2000,01,01,0,0,0); % starting date
pp = 1;
for i = 1:1096, % delta H for three years : size 1096 x 24
for j = 1:24,
    [y,m,d] = jd2date(st1+i-1);
    mjd(pp) = tmjd(y,m,d,j-1,30,0,1); %j-1 becasue kp strts from 0 Uhr 30 min - center of hourly mean
    F107(pp) = a11(i,2);
    kp(pp)=get_kp_new(mjd(pp)); 
    pp = pp+1;
end;
end;
%}